Nearest-Neighbor Distributions and Tunneling Splittings in Interacting Many-Body 

Two-Level Boson Systems 



O 

(N 

< 

(N 



Saul Hernandez-Quiro^l 

Instituto de Ciencias Fisicas, Universidad Nacional Autonoma de Mexico (UNAM), 62210 Cuernavaca, Morelos, Mexico and 
Facultad de Ciencias, Universidad Autonoma del Estado de Morelos (UAEM), 62209-Cuernavaca, Morelos, Mexico 

Luis BenetQ 

Instituto de Ciencias Fisicas, Universidad Nacional Autonoma de Mexico (UNAM), 62210 Cuernavaca, Morelos, Mexico 

(Dated: April 13, 2010) 

We study the nearest-neighbor distributions of the fc-body embedded ensembles of random matri- 
ces for n bosons distributed over two-degenerate single-particle states. This ensemble, as a function 
of k, displays a transition from harmonic oscillator behavior (k = 1) to random matrix type behavior 
(k = n). We show that a large and robust quasi-degeneracy is present for a wide interval of values of 
k when the ensemble is time-reversal invariant. These quasi-degenerate levels are Shnirelman dou- 
blets which appear due to the integrability and time-reversal invariance of the underlying classical 
systems. We present results related to the frequency in the spectrum of these degenerate levels in 
terms of k, and discuss the statistical properties of the splittings of these doublets. 
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I. INTRODUCTION 

The theoretical and experimental understanding of in- 
teracting many-body quantum systems has undergone 
considerable development in recent years. First, random 
matrix theory (RMT) has been quite successful in de- 
scribing the statistical properties of the fluctuations of 
the spectra of complex quantum systems, which include 
many-body interacting systems. Examples range from 
nuclear physics to disordered systems, including clasto- 
mechanical vibrations and quantum analog systems to 
classical chaotic billiards (see 0] for a detailed review). 
While this modeling has been quite successful, RMT is 
not a realistic theory since it assumes many-body forces 
between the constituents. More realistic stochastic model 
considering k body interactions are the embedded ensem- 
bles, initially introduced by Mon and French (2L This 
model can be defined for fermions and bosons 3|, and 
may be viewed as the generic models for stochasticity in 
many-body systems. 

Second, ultra-cold bosonic gases confined in optical lat- 
tices have become quite important due to the relative! 
simplicity to handle these systems experimentally 
In particular, Bose-Einstcin condensates (BECs) in a 
double- well potential is a common object of study Q. 
This system exhibits a great variety of interesting quan- 
tum phenomena, such as interference @, tunneling and 
self-trap ping 0, Q , Josephson oscillations [|| , and entan- 
glement (lOj . 

From the theoretical point of view, the two-level 
bosonic systems have been addressed using the mean field 
treatment of the Gross-Pitaevski equation 0[ll|, and the 



two-site Bose-Hubbard Hamiltonian. The latter can be 
written as 111, HI 



#bh = S(rn - m) - J(b\b 2 + blh) 
+ ^l n i( n i ~ 1) + n 2 (n2 - 1)] 



(1) 
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Here, b\ and hi are creation and annihilation operators for 
a boson on the ith site (i = 1, 2) and = S|Sj is the total 
number of bosons on that level, 5 is the energy difference 
of one-boson energies among the two sites, U is the on- 
site two-body interaction strength, and J is the hoping 
or tunneling parameter. The two-mode approximation 
in Eq. (JT|) is valid as long as the interaction energy U 
is much smaller than the level spacing of the external 
trap 0. 

The experimental observation of macroscopic tunnel- 
ing of bosons in a double well when the initial difference 
of population is below a critical value @, predicted in 
Ref. [7j, can be understood from the spectral properties 
of Eq. (|T|). For the simpler case S = 0, the spectrum 
consists of a lower region of nearly equidistant levels and 
an upper one displaying nearly degenerate doublets. The 
latter are actually responsible for the suppression of tun- 
neling; it has also been shown that coherences among 
nearby doublets yield oscillations with very small am- 
plitude [l4j . Taking the semiclassical limit, the system 
has a phase space representation similar to a pendulum, 
with the almost equidistant levels being associated with 
the libration zone and the nearly degenerate levels with 
the rotation zones. 

In this paper, we study the statistical properties of the 
spectrum of n bosons distributed on two levels coupled 
through random fc-body interactions. Thus, we merge 
the successful stochastic modeling of RMT with systems 
of the form of Eq. ([T]). This ensemble is actually a gen- 
eralization of the Bose-Hubbard type of Hamiltonians, 
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in particular with respect to the range of the interac- 
tion fc. Each member of the ensemble is Liouville in- 
tegrable (independently of fc) in the classical limit [Hj]. 
Yet, the spectral statistics of the ensemble correspond to 
a picket fence for k = 1, and follow RMT predictions for 
k = n [l(| . These facts make the ensemble somewhat spe- 
cial: completely integrable systems are associated with 
Poisson statistics, which is known as the Berry- Tabor 
conjecture [l7| . In addition, the spectral fluctuations 
of classically fully chaotic systems typically follow RMT 
predictions, which is known as the Bohigas-Giannoni- 
Schmit (or quantum-chaos) conjecture |l8(. We shall thus 
study the transition in the spectral statistics in terms of 
fc, considering the nearest-neighbor distribution as well 
as the occurrence and statistics of tunneling splittings. 
We shall address this for the case when the ensemble is 
time-reversal invariant (/3 = 1) or when this symmetry 
is broken (/? = 2). We find a systematic appearance of 
quasi-degeneracies on a large interval of k for the time- 
reversal case, which points out the underlying integrabil- 
ity properties of the members of the ensemble due to a 
theorem by Shnirelman (l9l [20j. Moreover, the number 
of such doublets as well as the statistics of the associated 
splittings display a dependence upon k. These results 
may be interesting for the understanding and modeling 
of three-body interactions in cold gases |2l[ . 

The paper is organized as follows. In Sec. UH we 
present the fc-body embedded ensembles of random ma- 
trices for two- level boson systems, and review some im- 
portant properties of this ensemble. In Sec. IIII1 we dis- 
cuss the nearest-neighbor distribution of the ensembles 
in terms of the interaction parameter k for both cases of 
Dyson's parameter f3. We obtain the systematic appear- 
ance of quasi-degenerate states in the spectrum linked to 
the /? = 1 case, and address the dependence of their num- 
ber with respect to k. In Sec. IIV1 we present the semi- 
classical limit of this ensemble and describe the structure 
of the corresponding phase space. Section [V] is devoted 
to the identification of the (3 = 1 quasi-degenerate states 
and present results on the statistical properties of their 
spacings. In Sec. lVIl we present a summary of our results 
and the conclusions. 



II. fc-BODY INTERACTING TWO-LEVEL 
BOSON ENSEMBLE 

We begin defining the most general fc-body interaction 
of n spin-less bosons distributed in two single-particle 
levels which, for simplicity, are assumed to be degenerate 
[case 5 = in Eq. (JTJ)] _ The single-particle states are asso- 
ciated with the operators bj and bj, with j = 1, 2, which, 
respectively, create or annihilate one boson on the single- 
particle level j. These operators satisfy the usual com- 
mutation relations for bosons. The normalized n-boson 
states are specified by = (M^ n) )- 1 (b\) r (bl) n - r \0), 

where Ay = [r\(n — 7*)!] 1 / 2 is a normalization constant 
and |0) is the vacuum state. The Hilbert-space dimen- 



sion is N = n + 1. In second-quantized form, the most 
general Hamiltonian Sfi with fc-body interactions can 
be written as 123 



6F = 



k 

E 

r.s—O 



(2) 



Physically, in Eq. ([2]) corresponds to n bosons con- 
fined, e.g., in a double — well potential, coupled only 
through fc-body interactions. Clearly, the degenerate 
Bose-Hubbard model Eq. ([T]) is a particular choice of the 

parameters for the combination H^ 1 2 1 + H 



k=2- 



Stochasticity is built into the Hamiltonian H k 



at the 

level of the fc-body matrix elements v;?J . These ma- 
trix elements are assumed to be Gaussian distributed 
independent random variables with zero mean and con- 

Then, v$v ( J?l, = v%(5 r >, s 5 r , s i + 



stant variance v% = 1. 



&p,iSr,r'&s,s')i where the over line denotes ensemble aver- 
age. As in the case of the canonical random matrix en- 
sembles [lj, Dyson's parameter /? distinguishes the sym- 
metry properties with respect to time-reversal invariance: 
P = 1 corresponds to the case where time-reversal sym- 
metry holds while the broken time-reversal case is de- 
noted by (3 = 2. The fc-body interaction matrix i;W 
is thus a member of the Gaussian orthogonal ensemble 
(GOE) for j3 = 1 or Gaussian unitary ensemble (GUE) 
for /3 = 2. This defines completely the fc-body embedded 
ensemble of random matrices for bosons distributed in 
I = 2 levels. Without loss of generality, in the follow- 
ing, we set i>o = 1. The combinatorial factors m 
Eq. are actually introduced in order to have an exact 
identity of the embedded ensembles with the canonical 
ensembles of RMT when fc = n 0, |22| . Indeed, the fac- 
tors cancel the square-root factors that appear by 
operating the fc = n creation and annihilation operators 
onto the many-body states l^i^}- Then, the central- 
limit theorem implies that the matrix elements become 
independent Gaussian distributed random variables; this 
is precisely the definition of the canonical ensembles of 
RMT. Consequently, for fc = n, all spectral fluctuations 
correspond exactly to the predictions of RMT. 

By construction, the number operator h — b\bi + b\i>2 

commutes with the Hamiltonian for all values of the 
rank of the interaction fc. The Hamiltonian is thus block 
diagonal in the occupation- number basis \^S n ') defined 
above. For a given value fc, the number of independent 
random variables of the ensemble is Kp{k) = /3(fc+l)(fc+ 
1 + ^3,i)/2, which in general is smaller than the Hilbert- 
space dimension N = n + 1. Therefore, for fc <§C n, the 
matrix elements of the Hamiltonian H^' are correlated, 
i.e., the number of independent matrix elements of the 
Hamiltonian is larger than the number of independent 
random variables. Moreover, some matrix elements are 
identically zero. 
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III. SPECTRAL STATISTICS IN TERMS OF k 

The evaluation of statistical measures of the spectrum 
requires unfolding the spectra, which removes the non- 
universal system-dependent contributions. This can be 
done by performing the unfolding individually for each 
spectrum (spectral unfolding) or by a single transfor- 
mation used for all members of the ensemble (ensemble 
unfolding). In the context of spectral fluctuations, er- 
godicity implies that the results arc independent of the 
unfolding method. 

In Ref. [22[, it was shown that the fc-body embedded 
ensemble of random matrices for bosons is non ergodic 
in the dense limit. The dense limit is defined as the limit 
n — > oo with k and the number of single-particle levels I 
fixed. This result was obtained analytically by consider- 
ing the fluctuations of the centroids and variances of in- 
dividual spectra, which do not vanish in the limit n — > oo 
of infinite Hilbert-space dimension Therefore, in the 
dense limit, ensemble average and spectral average yield 
in general different results. The non-ergodic character of 
the ensemble in the dense limit is a consequence of the 
fact that each member of the ensemble is Liouville inte- 
grable in the classical limit [15j . In this case, spectral 
unfolding is the only physically meaningful rectification 
of the spectra. In the numerical results described below, 
we implemented it by fitting the staircase function of each 
member of the ensemble separately with a polynomial of 
maximum degree 8. 

In Fig. [TJ we present the nearest-neighbor spacing dis- 
tribution Pfc(s) for various values of k, for (3=1 and 
n = 2000. These results were obtained after averaging 
over 1000 realizations of the ensemble. More details can 
be observed in the accompanying movie [H|. In these 
figures, we have included for comparison the Poisson dis- 
tribution and the Wigner surmise for the GOE 

For k = 1, the system corresponds to two coupled 
harmonic oscillators. Consequently, after unfolding, we 
obtain the expected distribution for an equidistant spec- 
trum, i.e., Pfe = i(s) = S(s — 1) [Fig.QJi]. As seen in Fig. [lb, 
for k = 2, this distribution changes considerably. It dis- 
plays a quite large peak at s = 0, a tail that decays 
somewhat slower than the Gaussian tail for larger values 
of s, and a broad peak around s = 1 reminiscent of the 
Dirac delta obtained for k = 1. The peak at s = in- 
dicates the occurrence of quasi-degenerate energy levels 
and, as we shall demonstrate below, it is a consequence 
of the time-reversal symmetry ((3 = 1) of the ensemble. 
Increasing slowly the value of k enhances the level clus- 
tering at s = and diminishes, shifts, and smoothes the 
peak at s = 1. This is illustrated for k = 10 in Fig. Q]:, 
where we also note that the tail of the distribution ap- 
proaches the exponential decay characteristic of the Pois- 
son distribution. The local maximum observed at s « 1 
disappears smoothly by increasing the value of k, being 
unnoticeable already for k = 75 [231 ]. 

By increasing the value of fc, the distribution Pk{s) 
evolves smoothly still displaying a strong degree of de- 



generacy at s — [cf. Figs. Id and le for k = 200 and 
k = 1000, respectively]. Eventually, around k = 1150, 
a new local maximum of the distribution is noticeable 
around s ~ 0.3 [Fig. [If], which moves toward larger val- 
ues of s for larger values of fc; this peak will become the 
single maximum of the GOE reached at k = n. From 
here on, except for the peak at s = 0, the distribution 
evolves toward the GOE results by increasing the value 
of k (see [23[), similarly to the transition observed in the 
spectral properties of the system when the dynamics of 
its classical analog evolves from near integrable to fully 
chaotic. Interestingly, the peak at s = is still observed 
for rather large values of k. Around k = 1850 [Fig. [Tg] 
this peak disappears, i.e., level repulsion completely sets 
in. Beyond k = 1900, the distribution corresponds essen- 
tially to that of a GOE. We note that there is no value 
of k where Pfc(s) fully coincides with the Poisson distri- 
bution, although it does so for larger spacings (tail of the 
distribution) in an extended range of values of k. It is not 
clear to us how to explain such an exponential tail for in- 
termediate values of k. At the moment, we believe that 
this fact may be related with a partial applicability of 
the original Berry- Tabor argument, which somehow can 
not be extended to all tori (see Ref. [24| for some recent 
results discussing the generic aspects of the Berry- Tabor 
conjecture). 

The remarkable property of the nearest-neighbor dis- 
tributions described above is the appearance and robust- 
ness of the large peak found around s = 0. This peak is 
not only pointing out the lack of level repulsion, but ac- 
tually indicating that a relevant part of the spectrum is 
degenerate or quasi-degenerate. This peak corresponds 
to the prediction of Shnirelman's theorem fl9j , which es- 
sentially states that smooth-enough time-reversal invari- 
ant ((3 = 1) and integrable Hamiltonian of two degrees of 
freedom have an asymptotically multiple spectrum, i.e., 
quasi-degenerate levels (see also dc}). Note that the as- 
sumptions of this theorem are fulfilled, since each mem- 
ber of the ensemble is Liouville integrable in the semi- 
classical limit [HI . 

To completely prove that the peak is indeed 
Shnirelman's peak, it suffices to consider the nearest- 
neighbor distribution Pk(s) for an ensemble of Hamiltoni- 
ans Hrf with broken time-reversal invariance, i.e., (3 = 2. 
If time- reversal is important, the peak should disappear 
for (3 = 2. The results are illustrated in Fig.[2]for different 
values of k, considering n = 1000 bosons and 1000 real- 
izations of the ensemble (see the corresponding movie [23| 
for more details). The figures show the transition from 
a picket fence spectrum (k = 1) to a GUE (k = n). In 
particular, they show the absence of the strong peak at 
,s = (Shnirelman's peak), even though as a function of 
k there is certain degree of level clustering, which are not 
quasi-degeneracies of the type discussed above. This is 
further illustrated in Fig. [3j where we show the relative 
number of levels fj,k corresponding to the first four bins 
of -Pfc(s) as a function of k, both for [3 = 1 and (3 = 2. 
In this figure we have also included the average number 
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FIG. 1: Nearest-neighbor distribution Pk{s) for the fc-body interacting two-level boson ensemble for /3 = 1, n = 2000 and 
(a) k = 1, (b) k = 2, (c) k = 10, (d) k = 200, (e) k = 1000, (f) k = 1150, (g) k = 1850 and (h) k = 2000. Notice the large peak 
observed at s = 0, which is linked with the occurrence of quasi-degenerate levels. The dashed curve corresponds to the Wigner 
surmise for j3 = 1, while the dotted curve is the Poisson distribution. 
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FIG. 2: Same as Fig. [T] for = 2, n = 1000 and (a) k = 1, (b) k = 2, (c) k = 10, (d) k = 200, (e) k = 500, (f) k = 700, 
(g) k = 800 and (h) k = 940. Notice that the strong peak observed in Fig. [T] at s = for /3 — 1 is absent in this case, indicating 
that its origin is due to time-reversal symmetry. Yet, certain degree of level clustering is still observed on a wide interval of k. 
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The classical Hamiltonian obtained in this way can 
be written as H% 1 2 , fa, fa) = Hok\h,h) + 

Vfc^(ii, h, <f>i, fa)- Here, W.o k (h,I 2 ) is a Hamiltonian 
that depends on the action variables only and is therefore 
integrable, and the perturbing term Vk (Ii,l2, fa, fa) 
carries all the dependence upon the angles. The first 
term is associated with all the diagonal contributions of 
Hjf^ , while the second one corresponds to all off-diagonal 
contributions. These terms are explicitly given by [16| 

n of = £ J^ Vs{h s)Vk - sil2 -^ k ~ s) ' (4) 



FIG. 3: (Color online) Relative measure fi k of the number 
of levels contained within the first four bins of the nearest- 
neighbor spacing distributions as a function of k/n. The blue 
curve (dotted curve with squares) corresponds to the time- 
reversal invariant case (/9 = 1) and the red curve (dashed 
curve with triangles) corresponds to the broken time-reversal 
case (/3 = 2). The continuous black curve (full circles) rep- 
resents the average number of Shnirelman doublets (/3 = 1) 
obtained using the symmetry properties of the cigenfunctions. 
The inset shows details for small values of k. 

of degenerate levels, which were identified using symmet- 
ric or antisymmetric combinations of the corresponding 
eigenfunctions (cf. Sect. IV A"f . Figure [3] implies that, as 
a function of k, there are different statistical properties 
of the degenerate levels. This in turn suggests the use of 
the quasi-degenerate levels, i.e., the tunneling splittings, 
as a possible measure to test fc-body interactions in such 
integrable systems. 

IV. SEMICLASSICAL LIMIT AND THE 
CLASSICAL PHASE SPACE 

A. Semiclassical limit 

Following Refs. [ID, [25|, we write an appropriate semi- 
classical limit for the algebraic Hamiltonian H k , which 
will allow us to identify systematically time-reversal 
related symmetric or anti-symmetric combinations of 

eigenfunctions. To this end, we symmetrize first Hjf 
with respect to the ordering of the creation and anni- 
hilation operators by exploiting the commutation rela- 
tions among the bosonic creation and annihilation oper- 
ators, typically in the form b\b s — (6f6 a + b s b\ — 8 r ^ s )/2 
(r,s = 1,2). Then, we use Heisenbcrg's semiclassical 
rules [H] 

b\ — -> 7 r 1/2 exp(iftv), b r — -> II' 2 cxp(-ifa), (3) 

where fa is an angle and I r is its canonically conjugated 
momentum. We emphasize the fact that considering the 
two-level case (I = 2) implies that the classical associated 
Hamiltonian has two degrees of freedom. 



V k =g 2jvW ) ™ s[{ s-t)(fa-fa)] 

[Vt{h-±,s) + Vt{h-\,t)} 

[V k - S {h -\,k-s) + V k - S {h - 3, ft - *)]• (5) 

In Eqs. (HJ) and ([5]), V t {I, s) are polynomials of degree t 
on the variable I defined as 

t 

V t (I,s) = H[I-(s-i)}, (6) 

i=l 

with s a numerical coefficient satisfying s > t > 0. We 
notice that the time-reversal symmetry properties of the 
ensemble are reflected in the matrix elements v f}. 

The classical Hamiltonian is therefore a general 
polynomial of degree k on the product of the actions 
with random coefficients, modulated by cosine functions 
whose argument is fa — fa. For ,8 = 1, the matrix el- 
ements v^~^ are real random numbers. Hence, time- 
reversal symmetry is manifested through the symmetry 
under reflection of both angles, i.e., fa — > —fa for both 
r = 1,2. In the case (3 = 2, the matrix «(' 9=2 ) is com- 
plex Hermitian, and the matrix elements can be written 

as w^ -2 ^ = ^ I exp[w r . s ], with the random phases 
satisfying v r<s = —v StT for Hcrmiticity. Therefore, the 
phases u rjj for j3 = 2 can be included into the cosine 
functions, manifestly breaking the invariance under si- 
multaneous reflections. 

From Eqs. (@| and (0, the angle variables appear in the 
Hamiltonian H k only through the combination fa — fa- 
in terms of the phase space geometry, this specific depen- 
dence corresponds to one single resonance, which implies 
the integrability of the classical Hamiltonian H k ^ ■ More 
explicitly, we perform a canonical transformation to new 
action and angle variables using the generating function 
W = Kfa + J {fa — fa), and obtain h = K - J, I 2 = J, 
X = fa j and ip = fa — fa . Substituting these expressions 
in Eqs. ([4]) and ([5]), the transformed Hamiltonian depends 
only on the angle ip. Since the angle \ docs not appear in 
the transformed Hamiltonian, its canonically conjugated 
action K is a conserved quantity, K = I\ + I 2 = n + 1, 
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with n the number of bosons. Therefore, besides the 
conservation of the energy (the Hamiltonian is time in- 
dependent), we have a second constant of motion, K. It 
is easy to show that the Poisson bracket between K and 
Wjj. is zero, thus implying that the Hamiltonian is (Li- 
ouville) integrable. For a fixed value of K, the reduced 
Hamiltonian (J, K, ip) is a time-independent one de- 
gree of freedom system with one parameter, which is al- 
ways integrable. In the language of symplectic geometry, 
the reduced Hamiltonian Hjf (J, K, ip) is identical to its 
normal form. In these variables, the population imbal- 
ance, which is a relevant quantity in the BECs context, 
is given by z = (h - I 2 )/(h + I 2 ) = 1 - 2J/K. 

We finish this discussion mentioning how to relate the 
action-angle variables of the reduced system (J and ip) to 
the actual coordinates and momenta q r and p r . r = 1,2, 
of the two sing le-particle modes. This is done by a lifting 
procedure [25j : Integration of the equations of motion of 
the reduced system yields J(t) and ip(t). Undoing the 
canonical transformation gives all I r (t) and (p r {t). Then 
we use the harmonic expressions 

jy 2 exp[T«0r] = (q r ± ip r )/V2, (7) 

which relate I r (t) and (f> r (i) to local coordinates and mo- 
menta, and thus yield the usual representation of the 
motion of each mode. 



B. Phase space structure in terms of k 

The essential features of the classical dynamics can 
be easily visualized in Poincare sections of the reduced 
Hamiltonian. Since the reduced Hamiltonian is a one- 
degree of freedom system, this representation corre- 
sponds to the level curves K, ip) = const. There- 
fore, the motion of a given initial condition follows the 
closed curve that includes the initial conditions. The ex- 
plicit appearance of the square-root factors in Eq. ([5]) 
makes the classical phase-space bounded, i.e., J <E [0,K] 
and ip G [— 7r, 7r], and therefore has the topology of a 
sphere. These properties are consistent with the inter- 
pretation of the two-mode Hamiltonian as a spin system. 

The phase-space structure of the unperturbed part 
U f\j,K) is trivial since J remains constant. There- 
fore, in a (Mercator) representation using the ip- J plane, 
the phase space appears foliated in horizontal straight 
lines for all fc, each representing a different level curve. 
The perturbing term vjf\j, K, ip) induces new structure 
due to the resonance; the argument of the cosine terms 
are fc-dependent integer multiples of ip. In Figs. [4) we 
present the phase-space structure for various values of k 
for (3 = 1 and n = 100; further details of the transition 
are given in the movie [23[ . While the system is integrable 
for all k and therefore the phase space is foliated by in- 
variant tori, the complexity of such tori increases with 
k. Figure S] was constructed choosing a specific (fixed) 
random matrix i/' 3 ' of dimension n+ 1, which defines the 



case k = n. For k = n — 1, we have defined the corre- 
sponding fc-body interaction by using the same matrix 
elements i>rf s for r, s < k, setting the remaining matrix 
elements to zero. This procedure can be iterated to ob- 
tain the corresponding matrix v with any desired value 
of k. 

For k = 1 [Fig. 0K| , we have two harmonic-oscillator 
wells centered around the two stable fixed points of the 
system, namely, around "0 = for small values of J and 
ip = ±7r at large values of J. The invariant curves as- 
sociated to initial conditions around such wells are self- 
retracing under time- reversal, i.e., each one is mapped 
onto itself under the transformation ip — > —ip. At inter- 
mediate values of J, the level curves are smooth deforma- 
tions of the unperturbed invariant curves (straight lines), 
thus illustrating Kolmogorov-Arnold-Moscr (KAM) the- 
orem. These tori are all self-retracing under time- 
reversal. 

For k = 2 [Fig. Hb] , the two harmonic fixed points (at 
the poles of the phase-space sphere) are still observed. 
Yet, at intermediate values of J, we notice two other fixed 
points close to ±7r/2 and the associated KAM-tori sur- 
rounding them. These structures are non-self-retracing 
under time reversal. Consequently, those tori that satisfy 
the Einstein-Brillouin-Keller (EBK) quantization rule 

S{E i ) = ±-jj{E i )A' l p = ±-{K i + ^), (8) 

where Ki is an integer and cti is the associated Maslov in- 
dex [27]], yield two degenerate ladders of levels, each pair 
associated with a torus and its associated time-reversed 
partner; quantum effects lift the degeneracy and produce 
the splitting of the quasi-degenerate levels. These dou- 
blets are precisely Shnirelman doublets. In Fig. Ht, we 
present the case k = 4, which displays the appearance of 
new non-self-retracing wells which yield ladders of double 
degenerate levels. Note that there are also self-retracing 
tori. The ladders of levels associated with these may dis- 
play accidental degeneracies with the levels of other lad- 
ders. These accidental degeneracies do not correspond 
to Shnirelman doublets since they may exist for the case 
of broken time-reversal invariance; this explains the level 
clustering observed for = 2 (see Fig.[3j). 

Increasing further the value of k increases the com- 
plexity of phase space [23| ■ New self- retracing and non- 
self-retracing wells appear in phase space, leading to 
Shnirelman doublets and perhaps accidental degenerate 
levels, respectively [see Fig. 3Jl]. Around k > 20, the 
number of stable fixed points begins to increase more 
rapidly, namely, quadratically with respect to k. Increas- 
ing k seemingly leads to a clustering in the sub-tropical 
region around the equator (J w K/2) of the majority of 
the stable fixed points (see Fig. 0f for k = 56). Even- 
tually, around k w 70, the stable fixed points begin to 
migrate to the polar regions of the phase-space sphere, 
with essentially all of them in that region for k > 81 [see 
Figs.Hg andHh]. 

As illustrated in Figs. [H the width of the spectrum is 
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FIG. 4: (Color online) Phase space representation (level curves) of the reduced Hamiltonian (J, K,tp) for /3 = 1 and 
n = 100. (a) k = 1, (b) fc = 2, (c) fc = 4, (d) fc = 13, (e) fc = 29, (f) k = 56, (g) fc = 81, (h) k = 100. The continuous 
lines represent three groups of 10 levels, taken from either edge and from the center of the spectrum. The color (grey) palette 
displays the overall energy scale of TL/f - (J, K,tj}), while the other one displays the full spectrum in the same scale. Note 
the clear appearance of ladders of levels associated with certain organizing centers in phase space. As k is increased, the level 
curves spread over all the available phase space, while the width of the spectrum decreases. 
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FIG. 5: Log-log plot showing the growth of the number of 
stable fixed points in terms of k for n = 100. The straight 
lines included have slope equal to 1 or 2. Around k > 20 the 
initial linear growth rate becomes quadratic. 



not constant with respect to k [221. For small k, the 
width is of the order of the energy scale spanned by 
(J, K, t/)) [see the spectral palette of Figs. Hk-EJi] ; this 
allows to visually identify different regions of the phase 
space where the eigenvalues are located. However, for 
larger values of k, the width is much smaller in compari- 
son to the full classical energy interval. In fact, for k = n 
[Fig. HJi], the width is proportional to n 1//2 [l|, while the 
full energy range is several orders of magnitude larger. 
Indeed, Eqs. Q and ([5]) involve homogeneous polynomi- 
als of degree k on J and J spans the interval [0,n + 1], 
thus spanning a huge energy interval. In turn, violent 
oscillations are due to the angular dependence, which is 
linear on k. Note that for k = n [Fig.@}i], the level curves 
associated with the energies of the spectrum spread es- 
sentially over all the phase space in such a way that it 
is not possible to distinguish one from another anymore, 
independently of their position in the spectrum. 

The growth rate of the number of stable fixed points 
is illustrated in a log-log plot in Fig. [5j The plot dis- 
plays an initial linear growth of the number of stable 
islands, which beyond k sa 20 becomes quadratic in k. 
This figure provides a solid base for the heuristic argu- 
ments of Ref. [l6|. There, it was argued that an expected 
quadratic growth in k of the number of stable fixed points 
(now shown in Fig. [5| diminishes the available phase- 
space area around the center of the wells. This implies 
that the EBK states which originally were found around 
those wells, will now be defined on tori which arc spread 
over more extended regions in phase space. Notice that 
this argument also shows that it is more difficult to have 
quantized states associated with non-self-retracing orbits 
for large values of k and thus explains that beyond cer- 
tain k, the number of Shnirclman pairs diminishes fast 
and eventually vanishes. 



V. SHNIRELMAN DOUBLETS AND THE 
STATISTICS OF THEIR SPLITTINGS 

As shown above, for /3 = 1 Slmirelman doublets appear 
and correspond to the quantization of non-self-retracing 
periodic orbits, i.e., orbits which are not mapped onto 
their selves under time- reversal invariance (ip — > —ip). 
Yet, these are not the only quasi-degenerate levels found 
since there are also accidental degeneracies involving two 
distinct self-retracing tori that just happen to have the 
same energy. Note that Shnirclman doublets disappear 
for f3 — 2, while the accidental degeneracies persist. 
Therefore, in order to distinguish the true Shnirclman 
doublets, we must consider the corresponding cigenfunc- 
tions, in a representation where the time-reversal in- 
variance is appropriately manifested. For this purpose, 
we first analyze the structure of the eigenfunctions of 
the quasi-degenerate states using a plane-wave decom- 
position which is straightforward to interpret in semi- 
classical terms [25j. Once we have classified the quasi- 



degeneracies, we address the question of the statistics of 
the Shnirelman splittings. 



A. Plane-wave decomposition of the wave functions 
and time-reversal symmetry 

The Hamiltonian Eq. ([3]) is conveniently expressed in 
the number occupation basis or Fock basis. For a given 
number of bosons n, we denote by |ni, ri2) the state hav- 
ing n\ bosons in the first single-particle state and n-i 
bosons in the second and n = ?il + n2. Upon diagonal- 
ization, the eigenfunctions of each realization of the en- 
semble are written as linear combinations of these basis 
states and have the form | $ r ) = J2 ni + n2 = n c ni,n2\ n i> n z)- 

The idea now is to use a representation where the sym- 
metry properties of the time-reversal invariance are mani- 
fested. To this end we recall that semiclassically the num- 
ber states can be represented as plane waves on the con- 
figuration torus (defined by the angle variables), namely, 
\ni) exp(zrii</>i )|</>i, </>2) [Hj]. Hence, the eigenstates 



can be written as |$ r 



+n 2 =n 



C ral,n2 



exp[i(ni0i 



ri 202)]|</ , i, </>2)- Since the total boson number is con- 
served, the associated dimensional reduction is imple- 
mented with the same canonical transformation de- 
scribed for the classical action-angle variables, which is a 
point transformation. Then, the eigenfunction of the rth 
excited state is written as 



$ r (V0 = exp(mx)^< 



%2 exp(in 2 t/j). (9) 



In Eq. ([5]), the factor exp(inx) is a common phase factor 
for all eigenstates, which can therefore be ignored, im- 
plying that the wave functions are functions only of the 
angle tp. Equation ([9]) defines the reduced representation 
of rth wave function. 

As discussed above, Shnirclman doublets are related 
to non-self-rctracing tori (under the transformation ip — ¥ 
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FIG. 6: (Color online) (a) Classical phase space represen- 
tation for k = 3 showing two pairs of non-self-retracing tori. 
The labels indicate the symmetry-related eigenfunctions. The 
color (grey) code is related to the full classical energy-scale 
of this case, (b) Reduced representation of the square mod- 
ulus of the eigenfunction A. (c) Linear combinations (|l(jp 
and with respect to levels A and A'; the figure shows 
that the linear combinations lie on opposite sides of the ip 
axis and therefore the eigenfunctions are a Shnirelman dou- 
blet, (d) Result of the linear combinations when considering 
two non-symmetry related nearby levels A and B. 



—ip) that sustain a state. The corresponding eigen- 
functions $r(?A) and $ r '(f/') appear as mixtures of wave 
functions localized around each symmetry-related torus. 
Therefore, we consider the linear combinations 

$+ r ,0/0 = -^Re($ r (V0+»MV0) ( 10 ) 

$- r ,(V0 = -^Re($ r (vo-iMVO). ( n ) 

In order to identify Shnirelman doublets we proceed 
as follows: First, we identify energy levels which lie very 
close together, which in practical terms means within the 
first few bins of the nearest-neighbor distribution mea- 
sured in units of the mean-level spacing. One would 
naively think that Shnirelman doublets appear as con- 
secutive levels; yet, accidental degeneracies due to other 
tori may have energies in between those of the doublets. 
This happens rather frequently for k/n > 0.2, where 
the number of stable fixed points grows quadratically on 
k. Therefore, we must check not only degeneracy with 
respect to the nearest level, but within a wider range. 
Then, for each candidate r of a Shnirelman doublet, we 
consider a second level r' and construct the linear com- 
binations given by Eqs. (fTU|) and ([TTj) . If and only if 
the functions $ J f r , (tp) and ^~ r ,(tp) are concentrated on 
one side of the ip axis (either positive or negative), and 
among them they are in opposite sides, then we say that 
the levels correspond to a Shnirelman doublet. In this 
case, the functions $* r ,(ip) and <&~ r ,{ip) are said to be 
related by the time-reversal transformation tp — > —ip. 

This method is illustrated in Fig. [5] In Fig. EJa), 
we plot the classical phase-space representation of the 
non-self-retracing tori corresponding to two pairs of 
Shnirelman doublets belonging to the same ladder. Fig- 
ures EJb)-|BJcL) display the reduced representation of the 
modulus square of some linear combinations involving 
these states. Figure displays one of the states cor- 
responding to the tori A. In Fig. Etc), we present the 
linear combinations (fTU| and (fTT]) involving the states A 
and A'; the resulting states are localized on either side of 
the ip = line, from where it is clear that these states are 
related by time reversal invariance. Finally, in Fig. EJd), 
we display the linear combinations involving two states 
A and B which belong to the same ladder but are not 
related by time-reversal invariance. 



B. Statistical properties of Shnirelman splittings 

We are interested in the spectral statistics of the spac- 
ings of the energies Ej and Ej> of Shnirelman doublets, 
i.e., the Shnirelman splittings, each one defines a charac- 
teristic tunneling time given by T w h/AEj. Let Ej and 
Ej' be the energies of a Shnirelman doublet of a specific 
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realization obtained by diagonalization. We define 



& = f , (13) 

where £' is the mean value of £J taken over the corre- 
sponding realization. Therefore, the quantity £ is a mea- 
sure of the splitting Ej and Eji in units of the energy 
of the doublet, averaged over the splittings of the cor- 
responding realization of the ensemble. The average is 
performed in order to compare the splittings of different 
realizations of the ensemble; hence, £ = 1. Then, £ is a 
kind of unfolded spacings of the Shnirclman doublets. 



In Fig. [71 we present the distribution of the spacings of 
Shnirelman doublets Pk{£) for various values of fc. These 
results were obtained for n = 2000 and 1000 realizations 
of the ensemble. For k = 2, this distribution displays a 
large narrow peak close to £ = 1.3xl0~ 9 which is not sym- 
metric. The peak decays very fast, acquiring a somewhat 
constant tail toward larger values of £; this tail vanishes 
after £ = 100. That is, while most doublets have a very 
small spacing in the normalized units used here, the dou- 
blets close to the edge of the ladder have larger spacings. 
For k = 3, similar results hold; yet, it is worth noticing 
that the distribution becomes somewhat wider with re- 
spect to the result for k = 2. As k increases further, the 
behavior of Pk (£) becomes more complex, with a gradual 
appearance of a second peak [close to £ = 1 in Fig. |7t] . 
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For a larger value of k around k « 750 [sec Fig. [7ji], both 
peaks have a similar amplitude; beyond this value of k, 
the left peak diminishes smoothly, eventually vanishing, 
yielding again a single peak distribution, this time for 
values centered around £ 1 . 

The transition in Pk (£) described above can be under- 
stood as follows. For small values or moderate values 
of k, the unimodal distribution reflects the existence of 
one or more ladders of Shnirelman doublets. Each lad- 
der has a number of doublets, the spacing with in each 
doublet becoming larger (smaller tunneling times) as we 
climb up the ladder. Imposing £ = 1 yields the long tail 
observed in the distribution. By increasing the value of 
k, the distribution Pfc(£) becomes bimodal, displaying a 
second peak centered around £ = 1. That is, the doublets 
have either a very small splitting or splittings of order 1. 
As mentioned above, larger splittings are attributed to 
the last doublets of a ladder around a stable fixed point. 
In order to have a significant number of them without 
increasing the number of small splittings, we conclude 
that their ladder must consist of very few (one or two) 
Shnirelman doublets. This idea is consistent with the 
fact that, for large value of k, only a single doublet is 
observed around the stable fixed points where the quan- 
tization condition Eq. ([5} holds. Note that the latter 
case implies again a unimodal distribution Pk(£), this 
time the peak being centered around 1, as it is observed 
numerically. 



VI. SUMMARY AND CONCLUSIONS 

In this paper, we have investigated the nearest- 
neighbor spacing distribution of the /c-body embedded 
ensembles for bosons distributed in two levels for /3 = 1 
and j3 = 2. For /3 = 1, we found a large peak at s = 
in a large interval of k which indicates the presence of 
degeneracies. This peak is quite robust in terms of k, 
disappearing only when k is very close to n, the total 
number of bosons. For (3 = 2, the peak is absent, despite 
of the fact that there are accidental quasi-degeneracies 
which yield small spacings; hence, the large peak is a con- 
sequence of the time-reversal invariance of the ensemble. 
We showed that this peak is associated with Shnirelman 
doublets, which semiclassically correspond to the quan- 
tization of tori that are non-self-retracing under time- 



reversal. These results provide further evidence on the 
integrability of the ensemble [IH based now on the spec- 
tral properties of the ensemble and therefore explain the 
non-ergodic properties of the ensemble [22[ • The fact that 
Shnirelman doublets are not observed for k very close to 
n, where GOE spectral statistics hold, is due to the fact 
that the non-self-retracing tori which would yield such 
doublets have a extremely small action, as shown in the 
phase-space representation of this case. 

We also found for /3 = 1 that the number of Shnirelman 
quasi-degeneracies displays a dependence upon k (cf. 
Fig. [3]). Moreover, the statistics of the normalized split- 
tings do display also a dependence on k; in particular, for 
k = 2 and k = 3 which are the physically relevant cases, 
we observe certain qualitative differences. We believe 
that these results may be interesting for understanding 
and modeling three-body interactions in Bose-Einstein 
condensates. 

Indeed, the existence of Shnirelman doublets opens 
the possibility of producing or observing other type 
of Josephson-like oscillations in two-mode Bose-Einstein 
condensates which may not be centered around zero pop- 
ulation imbalance [see, e.g., Figs.Sta) andHJb)]. To clar- 
ify this, we must emphasize that the degenerate states 
present in two-mode Bose-Einstein condensates 14| are 
not Shnirelman doublets; the degeneracies are due to the 
fact that the potential wells are indistinguishable, i.e., 
they are associated with the quantization of two tori re- 
lated by the symmetry J — > n + 1 — J. Therefore, in order 
to produce Shnirelman doublets, we must have the possi- 
bility of tuning all two-body interaction matrix elements 
at will, which may require considering also two species 
condensates. Once this is done, the statistical proper- 
ties of Shnirelman doublets could be used to characterize 
the role of interactions beyond k = 2. This will be the 
subject of a future work. 
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